index

1. Data

Show the code

# Library
library(dplyr)
library(ggplot2)

# load data 
## Isotope data 
isotope_data <-  utils::read.csv(here::here("data", "isotopic_data_2021.csv"), sep = ";", header = T,dec = ",")%>%
    # remove outlier values 
    filter(individual_code!="Arg-olf 42")

isotope_data_fish <- isotope_data %>%
  # only fish 
  filter(taxon == "Fish")%>%
  arrange(species)

## trawling data 
trawling_data_evhoe21 <-  utils::read.csv(here::here("data", "trawling_data_evhoe_2021.csv"), sep = ";", header = T, dec = ".")

# density distribution from biomass value 
density_distribution <- trawling_data_evhoe21%>%
  select(Nom_Scientifique, Tot_V_HV, Code_Station)%>%
  #selection of pelagic trawling 
  filter(Code_Station%in% c("Z0524", "Z0518", "Z0512", "Z0508", 
                            "Z0503", "Z0497", "Z0492"))%>%
  #selection of species sampled for isotope 
  filter(Nom_Scientifique%in%c("Arctozenus risso", "Argyropelecus olfersii", "Benthosema glaciale",
                               "Cyclothone", "Lampanyctus crocodilus", "Lampanyctus macdonaldi",
                               "Lestidiops sphyrenoides", "Maulisia argipalla", "Maurolicus muelleri",
                               "Melanostigma atlanticum", "Myctophum punctatum", "Notoscopelus bolini",
                               "Notoscopelus kroyeri", "Searsia koefoedi", "Serrivomer beanii",
                               "Xenodermichthys copei"))%>%
  mutate(Nom_Scientifique=recode(Nom_Scientifique, "Cyclothone"="Cyclothone spp."))%>%
  mutate(trawling_depth= case_when(Code_Station %in% c("Z0508") ~25,
                                   Code_Station %in% c("Z0492") ~370,
                                   Code_Station%in% c("Z0512") ~555,
                                   Code_Station %in% c("Z0503") ~715,
                                   Code_Station %in% c("Z0518") ~1000,
                                   Code_Station %in% c("Z0524") ~1010,
                                   Code_Station %in% c("Z0497") ~1335))%>%
  distinct()%>%
  group_by(Nom_Scientifique)%>%
  mutate(sum_sp=sum(Tot_V_HV))%>%
  ungroup()%>%
  group_by(trawling_depth, Nom_Scientifique)%>%
  mutate(pourcentage_bio=sum(Tot_V_HV/sum_sp*100))%>%
  # Selection of trawling depth 
  select(Nom_Scientifique, trawling_depth, pourcentage_bio)%>%
  # to have a round number to be able to multiply it afterwards 
  mutate(across(pourcentage_bio, round, 0)) %>%
  mutate(n_bio = as.integer(pourcentage_bio))%>%
  select(-c(pourcentage_bio))%>%
  tidyr::uncount(n_bio)

2. Isotopic niches

Ellipses

• ellipses at 40%
• Δ \(\delta\)¹³C = 2.36‰
• Δ \(\delta\)¹⁵N = 5.94‰

Show the code

niche_plot_community <- isotope_data%>%
  mutate(species= gsub("_"," ", species))%>%
  mutate(species=recode(species, "Cyclothone"="Cyclothone spp."))

# order by taxonomy 
niche_plot_community$species <- factor(niche_plot_community$species, 
                                       levels = c("Serrivomer beanii",
                                                  "Xenodermichthys copei",
                                                  "Maulisia argipalla",
                                                  "Searsia koefoedi",
                                                  "Cyclothone spp.",
                                                  "Argyropelecus olfersii",
                                                  "Maurolicus muelleri",
                                                  "Lestidiops sphyrenoides",
                                                  "Arctozenus risso",
                                                  "Benthosema glaciale",
                                                  "Lampanyctus crocodilus",
                                                  "Lampanyctus macdonaldi",
                                                  "Myctophum punctatum",
                                                  "Notoscopelus bolini",
                                                  "Notoscopelus kroyeri",
                                                  "Melanostigma atlanticum",
                                                  "Meganyctiphanes norvegica"))

# plot
ggplot(data = niche_plot_community, 
       aes(x = d13c, 
           y = d15n)) + 
  geom_point(aes(color = species, shape= taxon)) +
   paletteer:: scale_color_paletteer_d("pals::alphabet")+ 
   paletteer:: scale_fill_paletteer_d("pals::alphabet")+
  scale_shape_manual(values= c(19, 3))+
  scale_x_continuous(expression({delta}^13*C~'\u2030')) +
  scale_y_continuous(expression({delta}^15*N~'\u2030'))+
  stat_ellipse(aes(group = species, fill = species, color = species), 
               alpha = 0.2, level = 0.40,linewidth = 0.7, type = "norm", geom = "polygon")+
  theme_bw()+
  theme(legend.text = element_text(size=15),
        legend.title = element_text(size=15),
        axis.title = element_text(size=15),
        axis.text = element_text(size=15))+
  labs(shape="Taxon", col= "Species", fill="Species")+
  theme(aspect.ratio = 1)

Show the code

ggsave("niches_community.png", path = "figures", dpi = 700, height = 10, width = 12)

Overlaps

Show the code

# Prepare data 
mat_overlap_data <- isotope_data%>%
  filter(taxon=="Fish")%>%
  mutate(species= gsub("_"," ", species))%>%
  mutate(species=recode(species, "Cyclothone"="Cyclothone spp."))

# Arrange species by their d15n values 
mat_overlap <- mat_overlap_data%>%
  group_by(species)%>%
  mutate(mean_d15n=mean(d15n))%>%
  arrange(mean_d15n)%>%
  as.data.frame()

# ellipse at 40%
test.elp_community <- rKIN::estEllipse(data= mat_overlap,  x="d13c", y="d15n", group="species", levels=40, smallSamp = TRUE)

[1] "data.frame"
[1] "data.frame"
[1] "data.frame"
[1] "data.frame"
[1] "data.frame"
[1] "data.frame"
[1] "data.frame"
[1] "data.frame"
[1] "data.frame"
[1] "data.frame"
[1] "data.frame"
[1] "data.frame"
[1] "data.frame"
[1] "data.frame"
[1] "data.frame"
[1] "data.frame"

Show the code

# Extract the area of each polygon
elp.area_community <- rKIN::getArea(test.elp_community)

# determine polygon overlap for all polygons
elp.olp_community <- rKIN::calcOverlap(test.elp_community)%>%
  #renames row names
  mutate(OverlapID= gsub("_"," ", OverlapID))%>%
  mutate(OverlapID= gsub("40","", OverlapID))%>%
  mutate(across(where(is.numeric), round, 2))%>%
  tibble::column_to_rownames("OverlapID")%>%
  as.data.frame()

#renames colnames
colnames(elp.olp_community)<- gsub(x =colnames(elp.olp_community), pattern = "_", replacement = " ")
colnames(elp.olp_community)<- gsub(x =colnames(elp.olp_community), pattern = "40", replacement = "")

#plot 
ggcorrplot::ggcorrplot(elp.olp_community, lab = T, outline.color = "white", lab_size = 3, tl.cex = 10)+
  scale_fill_gradient2(limit = c(0,1), low = "white", high = "grey50", mid = "grey80", midpoint = 0.5)+
  labs(fill="Overlap value")+
  theme(axis.text = element_text(face="italic", size = 13),
        legend.text = element_text(size=10),
        legend.title = element_text(size=10),
        plot.background = element_rect(colour = "white"))

Show the code

ggsave("matrix_overlap.png", path = "figures", dpi = 700, width = 8, height = 6)

3. Depth segregation

Cluster

• input data for cluster = overlap matrix
• depth distribution = from complete 2021 trawling data (not only individuals sampled for isotope)

Show the code

# calculate median depth by species with trawling data set 2021
mean_depth_sp <- density_distribution %>%
  group_by(Nom_Scientifique) %>%
  mutate(median_depth = median(trawling_depth)) %>%
  select(Nom_Scientifique, median_depth) %>%
  distinct() %>%
  ungroup() %>%
  arrange(Nom_Scientifique) %>%
  rename(species = Nom_Scientifique)

# Gap statistic
factoextra::fviz_nbclust(
  elp.olp_community ,
  kmeans,
  nstart = 25,
  method = "gap_stat",
  nboot = 100,
  verbose = FALSE) +
  labs(subtitle = "Gap statistic method")

Show the code

res.km <- kmeans(scale(elp.olp_community), 5, nstart = 25)

# Dimension reduction using PCA
res.pca <- prcomp(elp.olp_community [,-5],  scale = TRUE)

# Individual coordinates
ind.coord <- as.data.frame(factoextra::get_pca_ind(res.pca)$coord)

# Add the clusters obtained using the k-means algorithm
ind.coord$cluster <- factor(res.km$cluster)

# Percentage of variance explained by dimensions
eigenvalue <- round(factoextra::get_eigenvalue(res.pca), 1)
variance.percent <- eigenvalue$variance.percent

# Add species names and median depth
ind.coord_depth <- cbind(ind.coord, mean_depth_sp)
ind.coord_depth$median_depth <- as.factor(ind.coord_depth$median_depth)

Cluster and depth distribution

Show the code

# dendrogram ----
dend <- elp.olp_community %>%
  dist() %>%
  hclust() %>%
  as.dendrogram()

#png("figures/dendrogram.png", units="in", width=6, height=4, res=700)
par(mar = c(1, 1, 1, 10))
dend %>%
  dendextend::set("labels_col",
      value = c("#86BBBD", "#ECA72C", "#4D85A8", "#9BABE8", "#D35D4A"),
      k = 5) %>%
  dendextend::set("branches_k_color",
      value = c("#86BBBD", "#ECA72C", "#4D85A8", "#9BABE8", "#D35D4A"),
      k = 5) %>%
  dendextend::set("labels_cex", 0.8) %>%
  dendextend::set("branches_lty", 2) %>%
  dendextend::set("leaves_bg") %>%
  dendextend::set("leaves_pch", 19) %>%
  dendextend::set("leaves_col",
    c(
      "#86BBBD",
      "#86BBBD",
      "#86BBBD",
      "#86BBBD",
      "#86BBBD",
      "#ECA72C",
      "#ECA72C",
      "#ECA72C",
      "#ECA72C",
      "#4D85A8",
      "#4D85A8",
      "#4D85A8",
      "#9BABE8",
      "#9BABE8",
      "#D35D4A",
      "#D35D4A"
    )
  ) %>%
  plot(horiz = TRUE, axes = FALSE)

Show the code

#dev.off()

# Density plot----
# assign each species to a cluster 
density_distribution_cluster <- density_distribution %>%
  mutate(
    cluster = case_when(
      Nom_Scientifique %in% c(
        "Argyropelecus olfersii",
        "Lampanyctus crocodilus",
        "Benthosema glaciale",
        "Myctophum punctatum",
        "Serrivomer beanii"
      ) ~ 4,
      Nom_Scientifique %in% c(
        "Lampanyctus macdonaldi",
        "Maulisia argipalla",
        "Searsia koefoedi"
      ) ~ 5,
      Nom_Scientifique %in% c(
        "Cyclothone spp.",
        "Notoscopelus bolini",
        "Notoscopelus kroyeri",
        "Melanostigma atlanticum"
      ) ~ 3,
      Nom_Scientifique %in% c("Xenodermichthys copei", "Maurolicus muelleri") ~
        1,
      Nom_Scientifique %in% c("Arctozenus risso", "Lestidiops sphyrenoides") ~
        2
    )
  ) %>%
  group_by(Nom_Scientifique) %>%
  arrange(desc(trawling_depth))

# Order in function of median depth
density_distribution_cluster$Nom_Scientifique = with(density_distribution_cluster, reorder(Nom_Scientifique, cluster, max))  

# plot
ggplot(density_distribution_cluster,
       aes(x = trawling_depth, y = Nom_Scientifique, group = Nom_Scientifique, 
           col=factor(cluster), fill=factor(cluster)))+ 
  scale_fill_manual(values = c("#D35D4A", "#9BABE8", "#ECA72C", "#86BBBD","#4D85A8"))+
  scale_color_manual(values = c("#D35D4A", "#9BABE8", "#ECA72C", "#86BBBD","#4D85A8"))+
  ggridges::stat_density_ridges(quantile_lines = TRUE, quantiles = 0.5 , alpha=0.4, size=0.7,
                                rel_min_height = 0.002, scale=1.2)+
  theme_bw()+
  scale_y_discrete(position = "left")+
  scale_x_reverse(limits = c( 1400,0))+
  coord_flip()+
  ylab(label = "")+ xlab("Depth (m)")+
  theme(axis.text.y = element_text(size=15),
        axis.text.x = element_text(face="italic", size=13, angle=80, vjust = 0.5, hjust=0.5),
        axis.title.x = element_text(size=15),
        axis.title.y = element_text(size=15))+
  guides(fill="none", col="none", alpha="none")

Show the code

ggsave("density_plot.png", path = "figures", dpi = 600, height = 8, width = 10)

4. Isotopic diversity index

Representativeness of sampling

• Percentage of species biomass in each station

Show the code

# Species % biomass and abundance by station
species_abundance <- trawling_data_evhoe21 %>%
  # selection of mesopelagic trawl
  filter(Code_Station %in% c("Z0524", "Z0518", "Z0512", "Z0508",
                             "Z0503", "Z0497", "Z0492")) %>%
  #deletion when only genus name and genus already sampled in isotopy + A. carbo
  filter(
    !Nom_Scientifique %in% c(
      "Cyclothone braueri",
      "Cyclothone microdon",
      "Myctophidae",
      "Aphanopus carbo"
    )
  ) %>%
  mutate(
    depth_layer = case_when(
      Code_Station %in% c("Z0508") ~ "epipelagic",
      Code_Station %in% c("Z0492", "Z0512") ~ "upper-mesopelagic",
      Code_Station %in% c("Z0503", "Z0518") ~ "lower-mesopelagic",
      Code_Station %in% c("Z0497") ~ "bathypelagic",
      Code_Station %in% c("Z0524") ~ "bottom-proximity"
    )
  ) %>%
  select(depth_layer, Nbr, Nom_Scientifique, Tot_V_HV) %>%
  group_by(Nom_Scientifique, depth_layer) %>%
  mutate(nb_ind = sum(Nbr)) %>%
  select(-Nbr) %>%
  distinct() %>%
  group_by(depth_layer) %>%
  mutate(sum_biomass_station = sum(Tot_V_HV)) %>%
  ungroup() %>%
  group_by(depth_layer, Nom_Scientifique) %>%
  mutate(pourcentage_biomass_sp = Tot_V_HV / sum_biomass_station * 100) %>%
  select(Nom_Scientifique, pourcentage_biomass_sp, nb_ind, depth_layer) %>%
  mutate(across(where(is.numeric), round, 2))
  
htmltools::tagList(DT::datatable(species_abundance))

Calculation

• Formatting of data and calculation of relative biomass within each station

Show the code

# Load data 
isotope_data_fish <- isotope_data %>% filter (species != "Meganyctiphanes_norvegica")

# sourcing the R functions from 'si_div' R script
source("R/si_div.R")

# Format indiviudal_si
individuals_si <- isotope_data_fish %>%
  select(individual_code, station, d13c, d15n, species_code) %>%
  rename(indiv_ID=individual_code,
         d13C= d13c,
         d15N =d15n,
         Species_code= species_code)%>%
  arrange(Species_code)

# species_status_biomass ----
# with all species sampled (not only species sampled for isotope)

# load total trawling data evhoe 2021
trawling_data_evhoe21 <-  utils::read.csv(here::here("data", "trawling_data_evhoe_2021.csv"), sep = ";", header = T, dec = ".")

species_status_biomass <- trawling_data_evhoe21 %>%
  # selection of mesopelagic trawl
  filter(Code_Station %in% c("Z0524", "Z0518", "Z0512", "Z0508",
                             "Z0503", "Z0497", "Z0492")) %>%
  #deletion when only genus name and genus already sampled in isotope + A. carbo
  filter(
    !Nom_Scientifique %in% c(
      "Cyclothone braueri",
      "Cyclothone microdon",
      "Myctophidae",
      "Aphanopus carbo",
      "Lampanyctus"
    )
  ) %>%
  # group by depth layer 
  mutate(
    Status = case_when(
      Code_Station %in% c("Z0508") ~ "epipelagic",
      Code_Station %in% c("Z0492", "Z0512") ~ "upper-mesopelagic",
      Code_Station %in% c("Z0503", "Z0518") ~ "lower-mesopelagic",
      Code_Station %in% c("Z0497") ~ "bathypelagic",
      Code_Station %in% c("Z0524") ~ "bottom-proximity"
    )
  ) %>%
  select(Status,
         Tot_V_HV,
         Nom_Scientifique,
         Code_Espece_Campagne) %>%
  distinct() %>%
  # sum species biomass by depth layer
  group_by(Nom_Scientifique, Status) %>%
  mutate(biomass_sp = sum(Tot_V_HV)) %>%
  select(-Tot_V_HV) %>%
  distinct() %>%
  # sum of total biomass by depth layer
  group_by(Status) %>%
  mutate(biomass_tot = sum(biomass_sp)) %>%
  # relative biomass of each species by station
  mutate(rel_biomass = biomass_sp / biomass_tot * 100) %>%
  select(-c(biomass_sp, biomass_tot)) %>%
  # selection of species sampled for isotopy in each depth
  filter(
    Status == "epipelagic" &
      Nom_Scientifique %in% c(
        "Lestidiops sphyrenoides",
        "Maurolicus muelleri",
        "Myctophum punctatum",
        "Notoscopelus kroyeri",
        "Notoscopelus bolini"
      ) |
      Status == "upper-mesopelagic" &
      Nom_Scientifique %in% c(
        "Arctozenus risso",
        "Argyropelecus olfersii",
        "Lampanyctus crocodilus",
        "Myctophum punctatum",
        "Notoscopelus kroyeri",
        "Xenodermichthys copei"
      ) |
      Status == "lower-mesopelagic" &
      Nom_Scientifique %in% c(
        "Arctozenus risso",
        "Argyropelecus olfersii",
        "Cyclothone",
        "Benthosema glaciale",
        "Lampanyctus crocodilus",
        "Maulisia argipalla",
        "Searsia koefoedi",
        "Serrivomer beanii",
        "Xenodermichthys copei"
      ) |
      Status == "bathypelagic" &
      Nom_Scientifique %in% c(
        "Argyropelecus olfersii",
        "Lampanyctus crocodilus",
        "Lampanyctus macdonaldi",
        "Myctophum punctatum",
        "Serrivomer beanii",
        "Xenodermichthys copei"
      ) |
      Status == "bottom-proximity" &
      Nom_Scientifique %in% c(
        "Argyropelecus olfersii",
        "Lampanyctus crocodilus",
        "Melanostigma atlanticum",
        "Serrivomer beanii",
        "Xenodermichthys copei"
      )
  ) %>%
  mutate(Species_code = tolower(Code_Espece_Campagne)) %>%
  rename(Species_name = Nom_Scientifique) %>%
  select(-Code_Espece_Campagne) %>%
  relocate(Status, .after = rel_biomass) %>%
  relocate(Species_code, .after = Species_name) %>%
  arrange(Species_name) %>%
  distinct()

htmltools::tagList(DT::datatable(species_status_biomass))

• percentage of biomass sampled for isotopy in each depth layer

Show the code

biomass_sampled <- species_status_biomass%>%
  group_by(Status)%>%
  mutate(sum_biomass = sum(rel_biomass))%>%
  select(sum_biomass, Status)%>%
  distinct()

htmltools::tagList(DT::datatable(biomass_sampled))

Epipelagic

_ 25m (Z0508)

Show the code

# 25m Z0508 ----
individuals_si_epi <- individuals_si%>%
  filter(station =="Z0508")%>%
  select(-station)

status_biomass_epi <- species_status_biomass%>%
  filter(Status=="epipelagic")

# computing mean Stable Isotope values for each species
# "group" column identical to species_code to fit with input format of function meanSI_group
# no "weight" input as number of indivuals sampled per species did not mirror actual species biomass
individuals_si_epi<-data.frame(group=individuals_si_epi[,"Species_code"], individuals_si_epi)
mean_si_species_epi<-meanSI_group(individuals_si_epi)

# computing coefficent of variation within each species to assess intraspecific variability
cbind(CV_d13C=mean_si_species_epi[,"sd_d13C"]/mean_si_species_epi[,"d13C"], CV_d15N=mean_si_species_epi[,"sd_d15N"]/mean_si_species_epi[,"d15N"] )

               CV_d13C    CV_d15N
lesti-sph -0.010766438 0.03266115
maur-mue  -0.005224151 0.05249590
myct-pun  -0.009814716 0.04047466
noto-bol  -0.009178283 0.02755544
noto-kro  -0.012155251 0.02163464

Show the code

# -> intraspecific variability is overall low (<20%)

# checking that species codes are the same in the two tables
row.names(mean_si_species_epi)==status_biomass_epi[,"Species_code"] # OK

     Species_code
[1,]        FALSE
[2,]         TRUE
[3,]         TRUE
[4,]         TRUE
[5,]         TRUE

Show the code

# building a single dataframe with all data for computing isotopic diversity indices
data_fish_epi <-data.frame(mean_si_species_epi[,c("d13C","d15N", "sd_d13C","sd_d15N")], rel_biomass=status_biomass_epi[,"rel_biomass"], Status=status_biomass_epi[,"Status"], latin_name=status_biomass_epi[,"Species_name"])

# scaling mean stable isotopes values using function "scale_rge01"
data_fish_scl_epi<-scaleSI_range01(data_fish_epi)

# computing isotopic diversity of the whole fish assemblage using scaled isotopic values and species relative biomass
ID_scl_ab_epi<-IDiversity(cons=data_fish_scl_epi, weight=data_fish_scl_epi[,c("rel_biomass")], nm_plot="1_epipelagic")

# printing results
result <- as.data.frame(round(ID_scl_ab_epi,3)) 

epipelagic_d13C_d15N

Upper mesopelagic

• 370m (Z0492) & 555m (Z0512)

Show the code

individuals_si_upm <- individuals_si%>%
  filter(station%in%c("Z0492", "Z0512"))%>%
  select(-station)

status_biomass_upm <- species_status_biomass%>%
  filter(Status=="upper-mesopelagic")

# computing mean Stable Isotope values for each species
individuals_si_upm<-data.frame(group=individuals_si_upm[,"Species_code"], individuals_si_upm)
mean_si_species_upm<-meanSI_group(individuals_si_upm)

# computing coefficent of variation within each species to assess intraspecific variability
cbind(CV_d13C=mean_si_species_upm[,"sd_d13C"]/mean_si_species_upm[,"d13C"], CV_d15N=mean_si_species_upm[,"sd_d15N"]/mean_si_species_upm[,"d15N"] )

              CV_d13C    CV_d15N
arct-ris -0.010421638 0.03438256
argy-olf -0.008706125 0.04283576
lamp-cro -0.015534945 0.07982239
myct-pun -0.019927081 0.04074265
noto-kro -0.013115883 0.02239182
xeno-cop -0.012364479 0.07722547

Show the code

# building a single dataframe with all data for computing isotopic diversity indices
data_fish_upm <-data.frame(mean_si_species_upm[,c("d13C","d15N", "sd_d13C","sd_d15N")], rel_biomass=status_biomass_upm[,"rel_biomass"], Status=status_biomass_upm[,"Status"], latin_name=status_biomass_upm[,"Species_name"])

# scaling mean stable isotopes values using function "scale_rge01"
data_fish_scl_upm<-scaleSI_range01(data_fish_upm)

# computing isotopic diversity of the whole fish assemblage using scaled isotopic values and species relative biomass
ID_scl_ab_upm<-IDiversity(cons=data_fish_scl_upm, weight=data_fish_scl_upm[,c("rel_biomass")], nm_plot="2_upper_meso")

# printing results
result <- as.data.frame(round(ID_scl_ab_upm,3)) 

370m

Lower mesopelagic

• 715m (Z0503)
• 1000m (Z0518)

Show the code

individuals_si_lwm <- individuals_si%>%
  filter(station%in%c("Z0503", "Z0518"))%>%
  select(-station)

status_biomass_lwm<- species_status_biomass%>%
  filter(Status=="lower-mesopelagic")

# computing mean Stable Isotope values for each species
individuals_si_lwm<-data.frame(group=individuals_si_lwm[,"Species_code"], individuals_si_lwm)
mean_si_species_lwm<-meanSI_group(individuals_si_lwm)

# computing coefficent of variation within each species to assess intraspecific variability
cbind(CV_d13C=mean_si_species_lwm[,"sd_d13C"]/mean_si_species_lwm[,"d13C"], CV_d15N=mean_si_species_lwm[,"sd_d15N"]/mean_si_species_lwm[,"d15N"] )

              CV_d13C    CV_d15N
arct-ris -0.013383792 0.02784671
argy-olf -0.010245107 0.03092443
bent-gla -0.015512985 0.06444602
cycl-otz -0.009264821 0.04946692
lamp-cro -0.020344506 0.04819340
maul-arg -0.009604780 0.03191368
sear-koe -0.022852986 0.05388211
serr-bea -0.012134196 0.06183851
xeno-cop -0.012017717 0.06196185

Show the code

# building a single dataframe with all data for computing isotopic diversity indices
data_fish_lwm <-data.frame(mean_si_species_lwm[,c("d13C","d15N", "sd_d13C","sd_d15N")], rel_biomass=status_biomass_lwm[,"rel_biomass"], Status=status_biomass_lwm[,"Status"], latin_name=status_biomass_lwm[,"Species_name"])

# scaling mean stable isotopes values using function "scale_rge01"
data_fish_scl_lwm<-scaleSI_range01(data_fish_lwm)

# computing isotopic diversity of the whole fish assemblage using scaled isotopic values and species relative biomass
ID_scl_ab_lwm<-IDiversity(cons=data_fish_scl_lwm, weight=data_fish_scl_lwm[,c("rel_biomass")], nm_plot="3_lower-mesopelagic")

# printing results
result <- as.data.frame(round(ID_scl_ab_lwm,3)) 

715m

Bathypelagic

• 1335m (Z0497)

Show the code

individuals_si_bathy <- individuals_si%>%
  filter(station=="Z0497")%>%
  select(-station)

status_biomass_bathy <- species_status_biomass%>%
  filter(Status=="bathypelagic")

# computing mean Stable Isotope values for each species
# "group" column identical to species_code to fit with input format of function meanSI_group
# no "weight" input as number of indivuals sampled per species did not mirror actual species biomass
individuals_si_bathy<-data.frame(group=individuals_si_bathy[,"Species_code"], individuals_si_bathy)
mean_si_species_bathy<-meanSI_group(individuals_si_bathy)

# computing coefficent of variation within each species to assess intraspecific variability
cbind(CV_d13C=mean_si_species_bathy[,"sd_d13C"]/mean_si_species_bathy[,"d13C"], CV_d15N=mean_si_species_bathy[,"sd_d15N"]/mean_si_species_bathy[,"d15N"] )

              CV_d13C    CV_d15N
argy-olf -0.009889227 0.01043410
lamp-cro -0.018677334 0.05649516
lamp-mac -0.022058396 0.02734756
myct-pun -0.024242014 0.04136229
serr-bea -0.014132768 0.05883833
xeno-cop -0.012585074 0.05902294

Show the code

# building a single dataframe with all data for computing isotopic diversity indices
data_fish_bathy <-data.frame(mean_si_species_bathy[,c("d13C","d15N", "sd_d13C","sd_d15N")], rel_biomass=status_biomass_bathy[,"rel_biomass"], Status=status_biomass_bathy[,"Status"], latin_name=status_biomass_bathy[,"Species_name"])

# scaling mean stable isotopes values using function "scale_rge01"
data_fish_scl_bathy<-scaleSI_range01(data_fish_bathy)

# computing isotopic diversity of the whole fish assemblage using scaled isotopic values and species relative biomass
ID_scl_ab_bathy<-IDiversity(cons=data_fish_scl_bathy, weight=data_fish_scl_bathy[,c("rel_biomass")], nm_plot="4_bathypelagic")

# printing results
result <- as.data.frame(round(ID_scl_ab_bathy,3)) 

bathypelagic

Near bottom

• Z0524

Show the code

# 1010 Z0524 ----
individuals_si_nb <- individuals_si%>%
  filter(station=="Z0524")%>%
  select(-station)

status_biomass_nb <- species_status_biomass%>%
  filter(Status=="bottom-proximity")

# computing mean Stable Isotope values for each species
individuals_si_nb<-data.frame(group=individuals_si_nb[,"Species_code"], individuals_si_nb)
mean_si_species_nb<-meanSI_group(individuals_si_nb)

# computing coefficent of variation within each species to assess intraspecific variability
cbind(CV_d13C=mean_si_species_nb[,"sd_d13C"]/mean_si_species_nb[,"d13C"], CV_d15N=mean_si_species_nb[,"sd_d15N"]/mean_si_species_nb[,"d15N"] )

              CV_d13C    CV_d15N
argy-olf -0.016533711 0.03529939
lamp-cro -0.015886548 0.04471856
mela-atl -0.009854421 0.04078758
serr-bea -0.012700557 0.05654679
xeno-cop -0.009943149 0.04105076

Show the code

# building a single dataframe with all data for computing isotopic diversity indices
data_fish_nb <-data.frame(mean_si_species_nb[,c("d13C","d15N", "sd_d13C","sd_d15N")], rel_biomass=status_biomass_nb[,"rel_biomass"], Status=status_biomass_nb[,"Status"], latin_name=status_biomass_nb[,"Species_name"])

# scaling mean stable isotopes values using function "scale_rge01"
data_fish_scl_nb<-scaleSI_range01(data_fish_nb)

# computing isotopic diversity of the whole fish assemblage using scaled isotopic values and species relative biomass
ID_scl_ab_nb<-IDiversity(cons=data_fish_scl_nb, weight=data_fish_scl_nb[,c("rel_biomass")], nm_plot="5_near_bottom")

# printing results
result <- as.data.frame(round(ID_scl_ab_nb,3)) 

bathypelagic

PCA

Show the code

# crate a data frame with all index value
trophic_indices <- data.frame(
  group = c(
    "Epipelagic",
    "Upper-mesopelagic",
    "Lower-Mesopelagic",
    "Bathypelagic",
    "Bottom proximity"
  ),
  IDiv = c(0.953, 0.862, 0.536, 0.932, 0.994),
  IDis = c(0.898, 0.823 , 0.491, 0.755, 0.318),
  IEve = c(0.448, 0.870 , 0.619, 0.633, 0.543),
  IUni = c(0.564, 0.800, 0.576, 0.713, 0.758)
)%>%
  tibble::column_to_rownames(var="group")

res.pca <- FactoMineR::PCA(trophic_indices, graph = FALSE)

factoextra::fviz_pca_biplot(res.pca, repel = TRUE,
                col.var = "#00778E", 
                col.ind = "gray50", 
                arrowsize = 1,
                title = ""
)

Show the code

ggsave("PCA.png", path = "figures", dpi = 700)

Show the code

htmltools::tagList(DT::datatable(trophic_indices))